#include "MCRadiationHeatFluxMaterial.h"
#include <iostream>
#include "FEProblem.h"
#include "libmesh/elem.h"
#include "libmesh/plane.h"
#include "libmesh/quadrature.h"
#include "Assembly.h"

using namespace std;
registerMooseObject("RoshanApp", MCRadiationHeatFluxMaterial);
template<>
InputParameters validParams<MCRadiationHeatFluxMaterial>()
{
  InputParameters params = validParams<Material>();
  params += validParams<RandomInterface>();
  params.addCoupledVar("temperature", "温度场");
  params.addRequiredParam<UserObjectName>("monte_carlo", "Monte Carlo");
  return params;
}

MCRadiationHeatFluxMaterial::MCRadiationHeatFluxMaterial(const InputParameters & parameters) :
      Material(parameters),
//	  RandomInterface(parameters, *parameters.get<FEProblem *>("_fe_problem"), parameters.get<THREAD_ID>("_tid"), false),
	_current_side_elem(_assembly.sideElem()),
	_temperature(coupledValue("temperature")),
	_MCuo(getUserObject<MCRadiationHeatFluxUserObject>("monte_carlo")),
	_radiation_heat_flux(declareProperty<Real>("MC_radiation_heat_flux")),
	_radiation_heat_flux_jacobi(declareProperty<Real>("MC_radiation_heat_flux_jacobi")),
	_recieve_radiation_heat_flux(declareProperty<Real>("MC_recieve_radiation_heat_flux")),
	_recieve_radiation_heat_flux_jacobi(declareProperty<Real>("MC_recieve_radiation_heat_flux_jacobi"))
{
}

void MCRadiationHeatFluxMaterial::computeQpProperties()
{
	Real sigma=5.67e-8;
	Real MC_flux = 0.0;
	Real MC_flux_jacobi = 0.0;
	Real Epsilon = 0.0;

	vector<Real> flux_radiation_list = _MCuo.getRadiationFluxList();
	vector<Real> flux_radiation_jacobi_list = _MCuo.getRadiationFluxJacobiList();
	vector<Point> center_point_list = _MCuo.getCenterPointList();
	vector<Real> epsilon_list = _MCuo.getEpsilonList();

	//for (unsigned int i = 0; i < flux_radiation_list.size(); ++i)
	//{
	//	cout << flux_radiation_list[i] << endl;
	//}

	int all_element_size = center_point_list.size();
	int findi = -1;

	for (unsigned int i = 0; i < all_element_size; ++i)
	{
		if ( ( center_point_list[i] - _current_side_elem->centroid() ).norm() < RaylineZero)
		{
			findi = i;
			break;
		}
	}

	if (findi == -1)
	{
		cout << "Current sideelement center : " << _current_side_elem->centroid() << endl;
		mooseError("Cannot find i.");
	}

	for(_qp = 0; _qp < _qrule->n_points(); ++_qp)
	{
		MC_flux = flux_radiation_list[findi];
		MC_flux_jacobi = flux_radiation_jacobi_list[findi];
		Epsilon = epsilon_list[findi];
		_radiation_heat_flux[_qp] = MC_flux - Epsilon * sigma * pow(_temperature[_qp], 4);
		_radiation_heat_flux_jacobi[_qp] = MC_flux_jacobi - 4 * Epsilon * sigma * pow(_temperature[_qp], 3);
		_recieve_radiation_heat_flux[_qp] = MC_flux;
		_recieve_radiation_heat_flux_jacobi[_qp] = MC_flux_jacobi;
	}
}
